- Vincent Fugère (il/lui)
- prof. en sciences de l’environnement à l’UQTR
- recherche en écologie des eaux douces au RIVE
- science des données & environnement
02-2022
Bonjour!
Aujourd’hui
- Analyses statistiques de base avec R
- Un peu de révision au passage:
dplyr,ggplot2 - Très peu de théorie! On assume que vous avez fait un cours universitaire de 1er cycle en statistiques et que vous êtes familier.e avec les termes “test de student”, “ANOVA”, “régression”…
- Masteriser la fonction
lm() - Code disponible ici
- Défis!
Données pour l’atelier
Nous allons réutiliser les ensembles de données que vous avez utilisé lors de deux atelier précédents: gapminder, palmer pinguins, starwars. Nous devons également charger les librairies dplyr et ggplot2.
library(dplyr)
library(ggplot2)
library(palmerpenguins)
library(gapminder)
data("starwars")
data("penguins")
data("gapminder")
# #au besoin
# install.packages(c("penguins","gapminder"))
Types de variables
- Selon StatCan
- variables catégoriques/qualitatives: nominales & ordinales
- variables numériques/quantitatives: continues & discrètes
head(penguins)
## # A tibble: 6 × 8 ## species island bill_length_mm bill_depth_mm flipper_length_… body_mass_g sex ## <fct> <fct> <dbl> <dbl> <int> <int> <fct> ## 1 Adelie Torge… 39.1 18.7 181 3750 male ## 2 Adelie Torge… 39.5 17.4 186 3800 fema… ## 3 Adelie Torge… 40.3 18 195 3250 fema… ## 4 Adelie Torge… NA NA NA NA <NA> ## 5 Adelie Torge… 36.7 19.3 193 3450 fema… ## 6 Adelie Torge… 39.3 20.6 190 3650 male ## # … with 1 more variable: year <int>
Défi!
Identifiez les variables qualitatives et quantitatives de l’ensemble de données starwars
Défi!
Identifiez les variables qualitatives et quantitatives de l’ensemble de données starwars
head(starwars)
glimpse(starwars)
## Rows: 87 ## Columns: 14 ## $ name <chr> "Luke Skywalker", "C-3PO", "R2-D2", "Darth Vader", "Leia Or… ## $ height <int> 172, 167, 96, 202, 150, 178, 165, 97, 183, 182, 188, 180, 2… ## $ mass <dbl> 77.0, 75.0, 32.0, 136.0, 49.0, 120.0, 75.0, 32.0, 84.0, 77.… ## $ hair_color <chr> "blond", NA, NA, "none", "brown", "brown, grey", "brown", N… ## $ skin_color <chr> "fair", "gold", "white, blue", "white", "light", "light", "… ## $ eye_color <chr> "blue", "yellow", "red", "yellow", "brown", "blue", "blue",… ## $ birth_year <dbl> 19.0, 112.0, 33.0, 41.9, 19.0, 52.0, 47.0, NA, 24.0, 57.0, … ## $ sex <chr> "male", "none", "none", "male", "female", "male", "female",… ## $ gender <chr> "masculine", "masculine", "masculine", "masculine", "femini… ## $ homeworld <chr> "Tatooine", "Tatooine", "Naboo", "Tatooine", "Alderaan", "T… ## $ species <chr> "Human", "Droid", "Droid", "Human", "Human", "Human", "Huma… ## $ films <list> <"The Empire Strikes Back", "Revenge of the Sith", "Return… ## $ vehicles <list> <"Snowspeeder", "Imperial Speeder Bike">, <>, <>, <>, "Imp… ## $ starships <list> <"X-wing", "Imperial shuttle">, <>, <>, "TIE Advanced x1",…
Le modèle linéaire général
- Une seule variable dépendante (réponse) quantitative
- Distribution normale (plus ou moins)
\[ y = mx + b \] \[ y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + ... + \beta_px_{ip} + \epsilon_i \] \[ \epsilon \sim N(0,\sigma^2) \]
Le modèle linéaire général
- Une seule variable dépendante (réponse) quantitative
- Distribution normale (plus ou moins)
hist(rnorm(10000))
Le modèle linéaire général
Test de Student (test de t)
- Différence entre deux moyennes
- Chez le manchot Adélie, est-ce que la masse corporelle diffère entre les mâles et les femelles?
penguins %>% filter(species == 'Adelie') %>% group_by(sex) %>% summarize(masse.moy = mean(body_mass_g))
## # A tibble: 3 × 2 ## sex masse.moy ## <fct> <dbl> ## 1 female 3369. ## 2 male 4043. ## 3 <NA> NA
Test de Student (test de t)
penguins %>% filter(species == 'Adelie', sex != 'NA') %>% ggplot(aes(x=body_mass_g, colour=sex)) + geom_density()
Test de Student (test de t)
Il semble y avoir une différence de masse entre les sexes. La fonction t.test nous permet de valider cette hypothèse.
adelie <- penguins %>% filter(species == 'Adelie', sex != 'NA') t.test(body_mass_g ~ sex, data = adelie)
## ## Welch Two Sample t-test ## ## data: body_mass_g by sex ## t = -13.126, df = 135.69, p-value < 2.2e-16 ## alternative hypothesis: true difference in means between group female and group male is not equal to 0 ## 95 percent confidence interval: ## -776.3012 -573.0139 ## sample estimates: ## mean in group female mean in group male ## 3368.836 4043.493
Test de Student (test de t)
Il semble y avoir une différence de masse entre les sexes. La fonction t.test nous permet de valider cette hypothèse.
adelie <- penguins %>% filter(species == 'Adelie', sex != 'NA') t.test(body_mass_g ~ sex, data = adelie)
## ## Welch Two Sample t-test ## ## data: body_mass_g by sex ## t = -13.126, df = 135.69, p-value < 2.2e-16 ## alternative hypothesis: true difference in means between group female and group male is not equal to 0 ## 95 percent confidence interval: ## -776.3012 -573.0139 ## sample estimates: ## mean in group female mean in group male ## 3368.836 4043.493
format(2.2e-16, scientific = F)
## [1] "0.00000000000000022"
Test de Student (test de t)
La fonction lm fonctionne pour tous les types de modèles linéaires (incluant le test de t).
modele <- lm(body_mass_g ~ sex, data = adelie) summary(modele)
## ## Call: ## lm(formula = body_mass_g ~ sex, data = adelie) ## ## Residuals: ## Min 1Q Median 3Q Max ## -718.49 -218.84 -18.84 225.00 731.51 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 3368.84 36.34 92.69 <2e-16 *** ## sexmale 674.66 51.40 13.13 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 310.5 on 144 degrees of freedom ## Multiple R-squared: 0.5447, Adjusted R-squared: 0.5416 ## F-statistic: 172.3 on 1 and 144 DF, p-value: < 2.2e-16
Test de Student (test de t)
\[ Y_i = \beta_0 + \beta_1x_{i} + \epsilon_i \] Quelles valeurs la variable “x” (sex) peut-elle prendre?
summary(modele)
## ## Call: ## lm(formula = body_mass_g ~ sex, data = adelie) ## ## Residuals: ## Min 1Q Median 3Q Max ## -718.49 -218.84 -18.84 225.00 731.51 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 3368.84 36.34 92.69 <2e-16 *** ## sexmale 674.66 51.40 13.13 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 310.5 on 144 degrees of freedom ## Multiple R-squared: 0.5447, Adjusted R-squared: 0.5416 ## F-statistic: 172.3 on 1 and 144 DF, p-value: < 2.2e-16
Défi!
Dans les données starwars, est-ce que les personnages provenant de Tatooine sont plus grands que ceux provenant de Naboo? On peut extraire les données pertinentes comme suit:
nab_tat <- starwars %>% filter(homeworld %in% c('Naboo','Tatooine'))
ggplot(aes(x=homeworld,y=height),data=nab_tat) + geom_boxplot()
Défi!
Dans les données starwars, est-ce que les personnages provenant de Tatooine sont plus grands que ceux provenant de Naboo?
t.test(height ~ homeworld, data = nab_tat)
## ## Welch Two Sample t-test ## ## data: height by homeworld ## t = 0.42386, df = 18.947, p-value = 0.6764 ## alternative hypothesis: true difference in means between group Naboo and group Tatooine is not equal to 0 ## 95 percent confidence interval: ## -22.27276 33.58185 ## sample estimates: ## mean in group Naboo mean in group Tatooine ## 175.4545 169.8000
Défi!
Dans les données starwars, est-ce que les personnages provenant de Tatooine sont plus grands que ceux provenant de Naboo?
modele <- lm(height ~ homeworld, data = nab_tat) summary(modele)
## ## Call: ## lm(formula = height ~ homeworld, data = nab_tat) ## ## Residuals: ## Min 1Q Median 3Q Max ## -79.455 -6.800 7.545 13.200 48.545 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 175.455 9.276 18.914 8.76e-14 *** ## homeworldTatooine -5.655 13.443 -0.421 0.679 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 30.77 on 19 degrees of freedom ## Multiple R-squared: 0.009227, Adjusted R-squared: -0.04292 ## F-statistic: 0.1769 on 1 and 19 DF, p-value: 0.6787
ANOVA
Lorsque notre variable qualitative indépendante comprend plus de deux niveaux, nous pouvons utiliser une analyse de variance (ANOVA) plutôt qu’un test de t. On ajuste le modèle de la même façon avec la fonction lm:
modele <- lm(body_mass_g ~ species, data = penguins) summary(modele)
## ## Call: ## lm(formula = body_mass_g ~ species, data = penguins) ## ## Residuals: ## Min 1Q Median 3Q Max ## -1126.02 -333.09 -33.09 316.91 1223.98 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 3700.66 37.62 98.37 <2e-16 *** ## speciesChinstrap 32.43 67.51 0.48 0.631 ## speciesGentoo 1375.35 56.15 24.50 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 462.3 on 339 degrees of freedom ## (2 observations deleted due to missingness) ## Multiple R-squared: 0.6697, Adjusted R-squared: 0.6677 ## F-statistic: 343.6 on 2 and 339 DF, p-value: < 2.2e-16
ANOVA
Les coefficients (béta) du modèle indiquent la différence de masse entre chaque niveau et le niveau de référence. Par défaut, le niveau de référence d’un facteur est celui qui arrive en première position lorsqu’on classe les niveaux en ordre alphabétique. Pour nos manchots, il y a trois espèces: Adelie, Chinstrap, Gentoo. Adelie est donc le niveau de référence.
levels(penguins$species)
## [1] "Adelie" "Chinstrap" "Gentoo"
coef(modele)
## (Intercept) speciesChinstrap speciesGentoo ## 3700.66225 32.42598 1375.35401
masse moyenne Adelie = 3700, Chinstrap = 3700+32, Gentoo = 3700+1375
ANOVA
masse moyenne Adelie = 3700, Chinstrap = 3700+32, Gentoo = 3700+1375
ggplot(aes(x=species,y=body_mass_g),data=penguins) + geom_boxplot()
ANOVA
Souvent, pour une ANOVA, on cherche plutôt à tester si la variable indépendante influence au moins un des niveaux, peu importe lequel, et sans nécessairement comparer les niveaux à un niveau de référence. On veut un tableau d’ANOVA qui ressemble à cela:
ANOVA
Il suffit d’utiliser la fonction anova pour obtenir le tableau d’ANOVA de notre modèle:
anova(modele)
## Analysis of Variance Table ## ## Response: body_mass_g ## Df Sum Sq Mean Sq F value Pr(>F) ## species 2 146864214 73432107 343.63 < 2.2e-16 *** ## Residuals 339 72443483 213698 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ANOVA
On complémente souvent le test de F de l’ANOVA avec des comparaisons par paires, pour déterminer quels groupes diffèrent les uns des autres. Dans ce cas, y a-t-il une différence entre Chinstrap et Adelie?
ggplot(aes(x=species,y=body_mass_g),data=penguins) + geom_boxplot()
ANOVA
Comparaisons par paires à l’aide de tests de Tukey (“Tukey’s Honest Significant Difference”):
TukeyHSD(aov(modele))
## Tukey multiple comparisons of means ## 95% family-wise confidence level ## ## Fit: aov(formula = modele) ## ## $species ## diff lwr upr p adj ## Chinstrap-Adelie 32.42598 -126.5002 191.3522 0.8806666 ## Gentoo-Adelie 1375.35401 1243.1786 1507.5294 0.0000000 ## Gentoo-Chinstrap 1342.92802 1178.4810 1507.3750 0.0000000
Défi!
À votre tour! Utilisez les données gapminder pour ajuster un modèle d’ANOVA. En 2007, y avait-il une différence d’espérance de vie (lifeExp) entre les continents (continent)? Le cas échéant, quelle(s) paire(s) de continents montraient une différence pour cette variable?
donnees_2007 <- gapminder %>% filter(year == 2007) ggplot(aes(x=continent,y=lifeExp),data=donnees_2007) + geom_boxplot()
Défi!
À votre tour! Utilisez les données gapminder pour ajuster un modèle d’ANOVA. En 2007, y avait-il une différence d’espérance de vie (lifeExp) entre les continents (continent)? Le cas échéant, quelle(s) paire(s) de continents montraient une différence pour cette variable?
modele <- lm(lifeExp ~ continent, donnees_2007) anova(modele)
## Analysis of Variance Table ## ## Response: lifeExp ## Df Sum Sq Mean Sq F value Pr(>F) ## continent 4 13060.7 3265.2 59.714 < 2.2e-16 *** ## Residuals 137 7491.2 54.7 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Défi!
À votre tour! Utilisez les données gapminder pour ajuster un modèle d’ANOVA. En 2007, y avait-il une différence d’espérance de vie (lifeExp) entre les continents (continent)? Le cas échéant, quelle(s) paire(s) de continents montraient une différence pour cette variable?
TukeyHSD(aov(modele))
## Tukey multiple comparisons of means ## 95% family-wise confidence level ## ## Fit: aov(formula = modele) ## ## $continent ## diff lwr upr p adj ## Americas-Africa 18.802082 13.827042 23.777121 0.0000000 ## Asia-Africa 15.922446 11.372842 20.472051 0.0000000 ## Europe-Africa 22.842562 18.155858 27.529265 0.0000000 ## Oceania-Africa 25.913462 11.183451 40.643472 0.0000305 ## Asia-Americas -2.879635 -8.299767 2.540497 0.5844901 ## Europe-Americas 4.040480 -1.495233 9.576193 0.2630707 ## Oceania-Americas 7.111380 -7.910342 22.133102 0.6862848 ## Europe-Asia 6.920115 1.763372 12.076858 0.0027416 ## Oceania-Asia 9.991015 -4.895221 24.877251 0.3464970 ## Oceania-Europe 3.070900 -11.857807 17.999607 0.9793776
Défi!
Merveilleux, mais… en regardant ces boîtes à moustaches (boxplot), trouvez-vous qu’il y a quelque chose de louche?
ggplot(aes(x=continent,y=lifeExp),data=donnees_2007) + geom_boxplot()
Suppositions des modèles linéaires
\[ Y_i = \beta_0 + \beta_1x_{i} + \epsilon_i \]
\[ \epsilon \sim N(0,\sigma^2) \]
Quelques suppositions des modèles linéaires:
- Les points de données sont indépendants (p. ex., on ne mesure pas le même manchot ou le même pays 10 fois!)
- Les résidus du modèle suivent une distribution normale
- La variance des résidus est identique pour tous les groupes (homoscédasticité)
Suppositions des modèles linéaires
Nous pouvons vérifier ces suppositions à l’aide de tests. IMPORTANT: pour ces deux tests, nous voulons que p > 0.05!
#1. normalité des résidus mod.resid <- resid(modele) # Sortir les résidus shapiro.test(mod.resid) # Test de Shapiro-Wilk: p < 0.05 = MAUVAIS!
## ## Shapiro-Wilk normality test ## ## data: mod.resid ## W = 0.96708, p-value = 0.001696
Suppositions des modèles linéaires
Nous pouvons vérifier ces suppositions à l’aide de tests. IMPORTANT: pour ces deux tests, nous voulons que p > 0.05!
#1. normalité des résidus mod.resid <- resid(modele) # Sortir les résidus shapiro.test(mod.resid) # Test de Shapiro-Wilk: p < 0.05 = MAUVAIS!
## ## Shapiro-Wilk normality test ## ## data: mod.resid ## W = 0.96708, p-value = 0.001696
# 2. homoscédasticité bartlett.test(lifeExp ~ continent, donnees_2007) # Test de Bartlett: p < 0.05 = MAUVAIS!
## ## Bartlett test of homogeneity of variances ## ## data: lifeExp by continent ## Bartlett's K-squared = 45.85, df = 4, p-value = 2.646e-09
Suppositions des modèles linéaires
Il faut vérifier les suppositions d’un modèle avant de l’interpréter! S’il y a un problème, on ne regarde même pas les résultats.
Solutions potentielles:
- exclure le(s) groupe(s) particulièrement problématique(s), ce qui diminue la portée de l’analyse
- transformer la variable réponse (p.ex. log)
- changer la structure du modèle en ajoutant des variables indépendantes
- utiliser un modèle linéaire généralisé ou d’autres outils statistiques plus avancés
ANOVA avec deux facteurs
Une ANOVA peut inclure plus d’une variables qualitatives, ce qui permet notamment de quantifier les intéractions entre deux ou plusieurs variables. Par exemple, est-ce que l’effet du sexe sur la masse des manchots dépend de l’espèce?
penguins <- tidyr::drop_na(penguins) ggplot(aes(x=species,y=body_mass_g,fill=sex),data=penguins) + geom_boxplot()
ANOVA avec deux facteurs
Une ANOVA peut inclure plus d’une variables qualitatives, ce qui permet notamment de quantifier les intéractions entre deux ou plusieurs variables. Par exemple, est-ce que l’effet du sexe sur la masse des manchots dépend de l’espèce?
modele <- lm(body_mass_g ~ species + sex, data = penguins) anova(modele)
## Analysis of Variance Table ## ## Response: body_mass_g ## Df Sum Sq Mean Sq F value Pr(>F) ## species 2 145190219 72595110 724.21 < 2.2e-16 *** ## sex 1 37090262 37090262 370.01 < 2.2e-16 *** ## Residuals 329 32979185 100241 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Que manque-t-il?
ANOVA avec deux facteurs
modele <- lm(body_mass_g ~ species * sex, data = penguins) anova(modele)
## Analysis of Variance Table ## ## Response: body_mass_g ## Df Sum Sq Mean Sq F value Pr(>F) ## species 2 145190219 72595110 758.358 < 2.2e-16 *** ## sex 1 37090262 37090262 387.460 < 2.2e-16 *** ## species:sex 2 1676557 838278 8.757 0.0001973 *** ## Residuals 327 31302628 95727 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ANOVA avec deux facteurs
N’oublions pas qu’il s’agit encore du modèle linéaire général:
\[ Y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + ... + \beta_px_{ip} + \epsilon_i \]
summary(modele)
## ## Call: ## lm(formula = body_mass_g ~ species * sex, data = penguins) ## ## Residuals: ## Min 1Q Median 3Q Max ## -827.21 -213.97 11.03 206.51 861.03 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 3368.84 36.21 93.030 < 2e-16 *** ## speciesChinstrap 158.37 64.24 2.465 0.01420 * ## speciesGentoo 1310.91 54.42 24.088 < 2e-16 *** ## sexmale 674.66 51.21 13.174 < 2e-16 *** ## speciesChinstrap:sexmale -262.89 90.85 -2.894 0.00406 ** ## speciesGentoo:sexmale 130.44 76.44 1.706 0.08886 . ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 309.4 on 327 degrees of freedom ## Multiple R-squared: 0.8546, Adjusted R-squared: 0.8524 ## F-statistic: 384.3 on 5 and 327 DF, p-value: < 2.2e-16
Régression & ANCOVA
Régression simple
Chez les personnages de star wars, y a-t-il une relation entre la masse corporelle et la grandeur?
ggplot(aes(x=height,y=mass),data=starwars) + geom_point()
Régression simple
Chez les personnages de star wars, y a-t-il une relation entre la masse corporelle et la grandeur?
ggplot(aes(x=height,y=mass),data=starwars) + geom_point()
Régression simple
Chez les personnages de star wars (autres que Jabba the Hut), y a-t-il une relation entre la masse corporelle et la grandeur?
sw <- starwars %>% filter(mass < 1000) ggplot(aes(x=height,y=mass),data=sw) + geom_point()
Régression simple
Chez les personnages de star wars (autres que Jabba the Hut), y a-t-il une relation entre la masse corporelle et la grandeur?
modele <- lm(mass ~ height, data = sw) summary(modele)
## ## Call: ## lm(formula = mass ~ height, data = sw) ## ## Residuals: ## Min 1Q Median 3Q Max ## -39.382 -8.212 0.211 3.846 57.327 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -32.54076 12.56053 -2.591 0.0122 * ## height 0.62136 0.07073 8.785 4.02e-12 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 19.14 on 56 degrees of freedom ## Multiple R-squared: 0.5795, Adjusted R-squared: 0.572 ## F-statistic: 77.18 on 1 and 56 DF, p-value: 4.018e-12
Régression simple
Validons les suppositions. Premièrement, la normalité des résidus:
mod.resid <- resid(modele) shapiro.test(mod.resid) # Test de Shapiro-Wilk: p < 0.05 = MAUVAIS!
## ## Shapiro-Wilk normality test ## ## data: mod.resid ## W = 0.912, p-value = 0.0004688
Ça ne passe pas…
Régression simple
On ré-essaie en transformant notre variable dépendante.
sw$log_masse <- log10(sw$mass) modele <- lm(log_masse ~ height, data = sw) mod.resid <- resid(modele) shapiro.test(mod.resid) #yé!
## ## Shapiro-Wilk normality test ## ## data: mod.resid ## W = 0.96702, p-value = 0.1156
Régression simple
Pour l’homogénéité de la variance, il suffit de regarder la distribution des résidus en fonction des valeurs de y estimées par le modèle.
Graphique tiré de ce site web
Régression simple
Pour l’homogénéité de la variance, il suffit de regarder la distribution des résidus en fonction des valeurs de y estimées par le modèle. Ici, on espère observer un nuage de point sans aucune forme triangulaire ou tendance directionnelle évidente.
plot(modele, which = 1)
Régression simple
Tout est bien. On peut interpréter le modèle…
summary(modele)
## ## Call: ## lm(formula = log_masse ~ height, data = sw) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.25043 -0.06977 0.01734 0.05502 0.22060 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.9911872 0.0690142 14.36 <2e-16 *** ## height 0.0048730 0.0003886 12.54 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.1052 on 56 degrees of freedom ## Multiple R-squared: 0.7374, Adjusted R-squared: 0.7327 ## F-statistic: 157.2 on 1 and 56 DF, p-value: < 2.2e-16
Régression simple
Tout est bien. On peut interpréter le modèle et le visualiser.
ggplot(aes(x=height,y=log_masse),data=sw) + geom_point() + geom_smooth(method='lm')
Régression simple
\[ y_i = \beta_0 + \beta_1x_{i} + \epsilon_i \]
Plutôt que d’utiliser ggplot pour ajuster le modèle à nouveau, nous pouvons aussi extraire les coéfficients (béta) de notre modèle et les utiliser pour “dessiner” la droite.
coef(modele)
## (Intercept) height ## 0.991187162 0.004872992
b0 <- coef(modele)[1] b1 <- coef(modele)[2] c(b0,b1)
## (Intercept) height ## 0.991187162 0.004872992
Quelle est l’équation de notre modèle linéaire?
Régression simple
\[ y_i = 0.9912 + 0.0049 * height_i + \epsilon_i \]
ggplot(aes(x=height,y=log_masse),data=sw) + geom_point() + geom_abline(intercept = b0, slope = b1)
Défi!
À votre tour! Dans les données penguins, y a-t-il une relation entre la longueur du bec (bill_length_mm) et la longueur des nageoires (flipper_length_mm)?
Défi!
modele <- lm(bill_length_mm ~ flipper_length_mm, penguins) par(mfrow=c(1,2)) plot(modele, which = 1:2)
Défi!
summary(modele)
## ## Call: ## lm(formula = bill_length_mm ~ flipper_length_mm, data = penguins) ## ## Residuals: ## Min 1Q Median 3Q Max ## -8.6367 -2.6981 -0.5788 2.0663 19.0953 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -7.21856 3.27175 -2.206 0.028 * ## flipper_length_mm 0.25482 0.01624 15.691 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 4.148 on 331 degrees of freedom ## Multiple R-squared: 0.4265, Adjusted R-squared: 0.4248 ## F-statistic: 246.2 on 1 and 331 DF, p-value: < 2.2e-16
longueur du bec = -7.21856 + 0.25482 * longueur de la nageoire
Défi!
ggplot(aes(x=flipper_length_mm,y=bill_length_mm),data=penguins) + geom_point() + geom_smooth(method='lm')
ANCOVA
Analyse de covariance: 1 variable qualitative, 1 variable quantitative.
ggplot(aes( x=body_mass_g, y=bill_length_mm, colour = species), data = penguins) + geom_point() + geom_smooth(method = "lm", aes(fill = species))
ANCOVA
ANCOVA
modele.sans.interaction <- lm(bill_length_mm ~ body_mass_g + species, data = penguins) summary(modele.sans.interaction) #equation pour le manchot Gentoo?
## ## Call: ## lm(formula = bill_length_mm ~ body_mass_g + species, data = penguins) ## ## Residuals: ## Min 1Q Median 3Q Max ## -6.8291 -1.6728 0.1244 1.5318 9.2904 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 24.908763 1.089730 22.858 < 2e-16 *** ## body_mass_g 0.003755 0.000289 12.991 < 2e-16 *** ## speciesChinstrap 9.908762 0.355289 27.889 < 2e-16 *** ## speciesGentoo 3.539179 0.499814 7.081 8.71e-12 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 2.419 on 329 degrees of freedom ## Multiple R-squared: 0.806, Adjusted R-squared: 0.8043 ## F-statistic: 455.8 on 3 and 329 DF, p-value: < 2.2e-16
ANCOVA
modele.avec.interaction <- lm(bill_length_mm ~ body_mass_g * species, data = penguins) summary(modele.avec.interaction) #equation pour le manchot Gentoo?
## ## Call: ## lm(formula = bill_length_mm ~ body_mass_g * species, data = penguins) ## ## Residuals: ## Min 1Q Median 3Q Max ## -6.4164 -1.6630 0.0683 1.5444 9.3138 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 27.1129056 1.6324193 16.609 < 2e-16 *** ## body_mass_g 0.0031599 0.0004371 7.228 3.48e-12 *** ## speciesChinstrap 5.0612873 3.3101846 1.529 0.127 ## speciesGentoo -0.5750259 2.7941191 -0.206 0.837 ## body_mass_g:speciesChinstrap 0.0013028 0.0008832 1.475 0.141 ## body_mass_g:speciesGentoo 0.0009698 0.0006225 1.558 0.120 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 2.414 on 327 degrees of freedom ## Multiple R-squared: 0.8081, Adjusted R-squared: 0.8051 ## F-statistic: 275.3 on 5 and 327 DF, p-value: < 2.2e-16
Régression multiple
Régression multiple
En ajoutant une troisième variable à notre modèle, on arrive à un modèle de régression multiple. Attention: les variables indépendantes ne doivent pas être trop corrélées les unes aux autres (problème de la multicolinéarité).
modele <- lm(bill_length_mm ~ body_mass_g + species + year, data = penguins) summary(modele)
## ## Call: ## lm(formula = bill_length_mm ~ body_mass_g + species + year, data = penguins) ## ## Residuals: ## Min 1Q Median 3Q Max ## -6.4997 -1.6727 0.0599 1.4895 9.5913 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -5.978e+02 3.270e+02 -1.828 0.0685 . ## body_mass_g 3.751e-03 2.879e-04 13.031 < 2e-16 *** ## speciesChinstrap 9.935e+00 3.541e-01 28.053 < 2e-16 *** ## speciesGentoo 3.540e+00 4.978e-01 7.110 7.28e-12 *** ## year 3.101e-01 1.629e-01 1.904 0.0578 . ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 2.41 on 328 degrees of freedom ## Multiple R-squared: 0.8082, Adjusted R-squared: 0.8058 ## F-statistic: 345.5 on 4 and 328 DF, p-value: < 2.2e-16
Régression multiple
On valide le modèle de la même façon.
par(mfrow=c(1,2)) plot(modele, which = 1:2)
Régression multiple
La représentation graphique d’un modèle de régression multiple est plus difficile. Souvent, nous nous contentons de montrer les coefficients du modèle sur un forest plot.
sjPlot::plot_model(modele)
Régression multiple
La représentation graphique d’un modèle de régression multiple est plus difficile. Souvent, nous nous contentons de montrer les coefficients du modèle sur un forest plot. Pour pouvoir comparer la taille des effets, il faut centrer-réduire toutes les variables quantitatives:
modele <- lm(bill_length_mm ~ scale(body_mass_g) + species + scale(year), data = penguins) sjPlot::plot_model(modele)
Au delà du modéle linéaire général
- La variable réponse ne suit pas une distribution normale: modèles linéaires généralisés
- Les points de données ne sont pas indépendants: modèles mixtes
- La relation entre x et y n’est pas linéaire: modèles additifs (GAM)
- Plusieurs variables dépendantes: analyse multidimensionnelle
Quelques ressources utiles: